;Given 16 function values f(i,j) for i,j=-1,0,1,2
;and normalized coordinates x(i),y(j) calculate
;coefficients of P(x)=SUM(i=0,3)SUM(j=0,3)a_ij*x^i*y^j
FUNCTION interp_bicube, f, x, y
  ;estimates derivative based on nearest two neighbors
  ; d_x( f(i,j) )~= ( f(i+1,j)-f(i-1,j) )/( x(i+1)-x(i-1) )
  fx_00=( f[2,1]-f[0,1] )/( x[2]-x[0] )
  fx_10=( f[3,1]-f[1,1] )/( x[3]-x[1] )
  fx_01=( f[2,2]-f[0,2] )/( x[2]-x[0] )
  fx_11=( f[3,2]-f[1,2] )/( x[3]-x[1] )
  ; d_y( f(i,j) )~= ( f(i,j+1)-f(i,j-1) )/( y(j+1)-y(j-1) )
  fy_00=( f[1,2]-f[1,0] )/( y[2]-y[0] )
  fy_10=( f[2,2]-f[2,0] )/( y[2]-y[0] )
  fy_01=( f[1,3]-f[1,1] )/( y[3]-y[1] )
  fy_11=( f[2,3]-f[2,1] )/( y[3]-y[1] )
  ; d_xy(f(i,j) )~= ( f(i+1,j+1)-f(i+1,j-1)-f(i-1,j+1)+f(i-1,j-1) ) / $
  ;                     ( ( x(i+1)-x(i-1) )( y(j+1)-y(j-1) ) )
  fxy_00=( f[2,2]-f[2,0]-f[0,2]+f[0,0] )/( (x[2]-x[0])*(y[2]-y[0]) )
  fxy_10=( f[3,2]-f[3,0]-f[1,2]+f[1,0] )/( (x[3]-x[1])*(y[2]-y[0]) )
  fxy_01=( f[2,3]-f[2,1]-f[0,3]+f[0,1] )/( (x[2]-x[0])*(y[3]-y[1]) )
  fxy_11=( f[3,3]-f[3,1]-f[1,3]+f[1,1] )/( (x[3]-x[1])*(y[3]-y[1]) )

  ;assemble the RHS vector
  f_vec=[ f[1,1], f[2,1], f[1,2], f[2,2],                                      $
           fx_00,  fx_10,  fx_01,  fx_11,                                      $
           fy_00,  fy_10,  fy_01,  fy_11,                                      $
          fxy_00, fxy_10, fxy_01, fxy_11 ] 

  ;assemble the inversion matrix (Wikipedia article on bicubic interpolation)
  ;       1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16
  A_1 = [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]  ;1
  A_2 = [ 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]  ;2
  A_3 = [-3, 3, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]  ;3
  A_4 = [ 2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]  ;4
  A_5 = [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]  ;5
  A_6 = [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]  ;6
  A_7 = [ 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0]  ;7
  A_8 = [ 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0]  ;8
  A_9 = [-3, 0, 3, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0]  ;9
  A_10= [ 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0,-2, 0,-1, 0]  ;10
  A_11= [ 9,-9,-9, 9, 6, 3,-6,-3, 6,-6, 3,-3, 4, 2, 2, 1]  ;11
  A_12= [-6, 6, 6,-6,-3,-3, 3, 3,-4, 4,-2, 2,-2,-2,-1,-1]  ;12
  A_13= [ 2, 0,-2, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0]  ;13
  A_14= [ 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 1, 0, 1, 0]  ;14
  A_15= [-6, 6, 6,-6,-4,-2, 4, 2,-3, 3,-3, 3,-2,-1,-2,-1]  ;15
  A_16= [ 4,-4,-4, 4, 2, 2,-2,-2, 2,-2, 2,-2, 1, 1, 1, 1]  ;16

  A_inv=[ [A_1],[A_2],[A_3],[A_4],[A_5],[A_6],[A_7],[A_8],                     $
          [A_9],[A_10],[A_11],[A_12],[A_13],[A_14],[A_15],[A_16] ]

  coeffs=TRANSPOSE(A_inv)#f_vec
  RETURN,coeffs
END
